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Abstract 

We describe the basic features of SELEN, an open source Fortran 90 program for the nu- 
merical solution of the so-called "Sea Level Equation" for a spherical, layered, non-rotating 
Earth with Maxwell viscoelastic rheology. The Sea Level Equation was introduced in the 
70s to model the sea level variations in response to the melting of late-Pleistocene ice- 
sheets, but it can be also employed for predictions of geodetic quantities such as vertical 
and horizontal surface displacements and gravity variations on a global and a regional scale. 
SELEN (acronym of SEa Level EquatioN solver) is particularly oriented to scientists at their 
first approach to the glacial isostatic adjustment problem and, according to our experience, 
it can be successfully used in teaching. The current release (2.9) considerably improves the 
previous versions of the code in terms of computational efficiency, portability and versatility. 
In this paper we describe the essentials of the theory behind the Sea Level Equation, the 
purposes of SELEN and its implementation, and we provide practical guidelines for the use of 
the program. Various examples showing how SELEN can be configured to solve geodynamical 
problems involving past and present sea level changes and current geodetic variations are also 
presented and discussed. 

Keywords: Sea Level Equation - Glacial Isostatic Adjustment - Relative Sea Level - 
Geodetic Variations 

1 Introduction 

The physical processes governing sea level changes of glacio-isostatic and hydro-isostatic ori- 
gin have been thoroughly studied in the last decades. They are described in the so-called 
"Sea Level Equation" (hence after referred to as SLE), the integral equation first introduced 
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bv lFarrell and Clarkl (Il976l) to model sea level variations following the melting of late-Pleistocene 
ice-sheets. The SLE is also currently employed to study the s ea level changes assoc iated with 
present terrestrial ice melting in response to global warming (e.g. Mitrovic a et al.ll200ll) , a critical 
issue for its potential impact on environment and society ( Solomon et all 20071) . The variation 
of currently observed geodetic quantities (3-D crustal deformations, gravity variations, altimetric 
and tide gauges signals) are directly affected by the melting of past and present ice sheets. How- 
ever, they are also sensitive to the meltwater component of the surface load, whose time evolution 
is governed by the same equation. The wide spectrum of possible applications of the SLE confirms 
the c entral role of glacial isostatic adjustment ( GIA) in the frame work of both satellite (jPeltierl 



2004 ) and ground-based geodetic investigations ( King et al. . 2010). 



SELEN is an open source program written in Fortran 90, primarily designed to simulate the 
sea level and geodetic variations in response to the melting of continental ice sheets. SELEN pro- 
vides a numerical solution to the SLE and gives access to various related quantities of interest 
for geodynamics and geo desy. The program, which was first introduced in a primitive form by 
Spada and Stocchl ( 20071) and now obsolete (hereafter referred to as SELEN 1.0), has recently 



been deeply re-organized and improved in a number of aspects, ranging from the structure of the 
program to the computational performance and portability. In particular, all components have 
been reviewed and optimized, and various utilities have been added to facilitate the numerical 
solution to problems of geodetic relevance, as the time-variations of the harmonic coefficients of 
the gravity field. Recently, SELEN 2.9 has been successfully employed, for practical applications 
in the context of GIA, within the 2nd Training School on Glacial Isostatic Adjustment Modeling, 
held in Gavle (Sweden) in June 2011. The feedback from the students and from a number of 
scientists involved in various fields of geodynamics and geodesy, who are rout inely working with 
SELE N, have greatly contributed to improve the code previously released by Spada and Stocchl 
(|2007j ) and has stimulated the release of this upgraded version. 

Though the principal physical ingredients of the SLE are implemented in SELEN, some ap- 
proximations are adopted. First, SELEN assumes a linear incompressible rheology and a spheri- 
cally symmetric undeformed Ea rth. Therefor e, the program ta kes advantage from the viscoelastic 



Green's function formalism (see lPeltierlll974l and references in lSpada et al.ll2011l ). Consequently, 



lateral rheological variations are not taken into account. As a second approximation, possible ef- 
fec ts of the Earth's rotat i onal fl uctuations upon sea level variations are not modeled. As discussed 
by iMilne and Mitrovica (Il998l) . the rotational effects are important in some particular regions 
of the globe, but not la rge enough to invalidate the whole solution of the SLE. Third, following 



Farrell and Clarkl (|1976I ). SELEN does not account for the horizontal migration of shorelines in 



response to sea level change. This indeed constitutes a crude approximation especially in areas of 
shallow bathymetry, which will be relaxed in the future rele ases of SELEN so to allow for "gravita- 
tionally self-consistent" paleogeographical reconstructions (Peltier, 2004) . Finally, it is important 
to remark that SELEN does not account for tectonic contributions to sea level, for local effects 
such as subsidence driven by sediment loading, nor any possible anthropic contribution to sea 
level change. Furthermore, SELEN does not account for ocean dynamics or possible steric sealevel 
variations, nor it accounts for ice dynamics. The challenging topic of the coupling between the 
SLE, the ocean circulation and the continental cryosphere is left to others. 

This manuscript is not self-contained. An exhaustive description of the technical details is 
avoided here, and emphasis is given on the general features of SELEN and the discussion of a 
set of case studies of genera l reley ance. The theory behind SELEN has been illustrated in detail 
bv lSpada and Stocchi ( 20061 2007 ). to which the reader is also referred for references to the most 
important papers published on this topic. SELEN does not come with an user guide. However, 
the comments across all Fortran program units and the self-explicatory SELEN configuration file 
should be sufficient. The paper is organized as follows. The physics behind SELEN is briefly 
summarized in Section [2] Section [3] illustrates how SELEN works. Section [4] introduces a few 
case studies where SELEN is employed to model past sea level variations (subsection 14. ip and 
present-day variations of geodetic quantities (|4.2[) . Finally, in Section [5] we draw our conclusions. 

SELEN is free software and anyone is welcome to distribute it under certain conditions within 
the terms of a Global Public License (GPL) (for details, visit http://www.gnu.org/licenses/). The 
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SELEN package is available upon request to the authors. 



2 Theory 



Here we briefly illustrate the background theory of the SLE, the integral equation that describes 
the sea level variations and solid Earth deformations associated with GI A. Essentially, the materia l 
that follows is a condensed summary of the SLE theory first exposed bv lFarrell and Clarkl (Il97fih . 
on which SELEN is based. In addition, some theory notes will also be provided in Section @] to 
discuss specific SELEN outputs. For furt her details, including t he numerical implementatio n of 
the so-called "pseudo-spectral" method (jMitrovica and Peltierl | l991f | M itrovic a et al. on 
which SELEN relies upon, the reader is referred to the review of Spada and Stocchi (|2006l) and to 
references therein. 

As it will appear from the synthesis below, the SLE (and its numerical implementation) does 
not involve absolute quantities and can, consequently, only provide variations of geophysical and 
geodetic quantities relative to a reference state. To clearly illustrate what SELEN is actually 
computing, it is convenient to first define "sea level" as 



SL(w,t) = R SS -R a 



(1) 



where u> = (8,X), 9 is colatitude and A is longitude, t is time, and with R ss and R se we denote the 
radius of the (equipotential) sea surface and of the solid surface of the Earth, both relative to the 
Earth's center of mass, respectively. Quantity SL(w,i), which has indeed an "absolute" character, 
is not what SELEN is solving for. Rather, the quantity involved in the SLE is sea level change 



(2) 



where t r is a reference time that, in the numerical applications of the SLE generally denotes 
a remote epoch prior to the beginning of ice melting. Hence, from ([T]), the definition of sea 
level change suitable for GIA studies is 



S(u,t) =N-U, 



(3) 



where N represents the sea surface variation and U is the vertical displacement of the solid surface 
of the Earth. Basically, Eq. Q represents the SLE in its simplest form; what follows is aimed to 
illustrate the relationship between S and the variations of the ice thickness through time, in order 
to obtain a form amenable to a numerica l approach. 

According to iFarrell and Clark (ll976lL N can be written as 



N(u,t) = G + c, 

where c is a yet undetermined function of time, and the geoid height variation is 

4> 



7 



(4) 



(5) 



in which 7 is the reference gravity at the surface of the Earth and $(w, t) is the total variation of 
the gravity potential. Hence, using Eq. {5} into © gives 



S(w,t) = U + c, 

7 

where mass conservation of the system (Ice sheets + Oceans) is ensured taking 

JTlj(t) 



*-Z7 

7 



(6) 



(7) 
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where the density of water p w is assumed to be constant, rrii is the mass variation of the ice sheets, 
A a is the (constant) area of the present-day oceans and the overbar indicates the average over the 
surface of the oceans 

T^j=±-J(...)dA, (8) 

where dA = a 2 sin 9d9dX is the area element and a is Earth average radius. From Eq. ([5]), the SLE 
can be therefore written as 



S{ U ,t)=(Z-U)+S*-(Z-Uy (9) 
where the "eustatic" sea level variation 

S B {t) = _n (1Q) 

shows the remarkable property S E = S. The SLE has solution S — S only in the case of a rigid, 
non self-gravitating Earth ([/ = $ = in Eq. [9]). 

Functions U (u, t) and t) will depend on the spatiotemporal variations of the surface load. 
This is expressed by 

C(u,t) = Pi I + p w SO, (11) 

where the two terms on the right hand side are associated with the waxing and waning of the ice 
sheets, and with the redistribution of meltwater in the ocean basins, respectively. In Eq. (jlll) . pi 
is ice density, O is the "ocean function" (O = 1 on the oceans, and O = on land), and 

I{u,t)=T-T , (12) 

is the ice thickness variation, where T(u,t) is absolute ice thickness, and Tq(uj) is a reference 
thickness (e.g. the thickness at the Last Glacial Maximum, LGM, 21 kyrs ago). The mass 
variation in Eq. (110[) is obtained from (|12[) by integration over the ice-covered regions: 



rm(t) = J pil dA. (13) 

According to Eq. (fTTj) . vertical displacement stems from two terms 

U(u,t) = piG u ®iI + p w G u ® S, (14) 

where G u is the Green's function for vertical displacements, ®i and ® are spatio-temporal con- 
volutions over the ice- and ocean-covered regions, respectively. Similarly 

= p l G 4> <Z> l I + puGptgio S, (15) 

where G^ is the Green's function for increm e ntal g ravitational potential. Explicit expressions for 
G u and are given in ISpada and Stocchi ( 2006 ) in terms of the load-deformation coefficients 



(LDCs) h{t) and k(t), respectively. Introducing the sea level Green's function G s as 

^l(u;,t) = ^-G u , (16) 

7 7 

substitution of Eqs. (fT5|) and (fTU) into © gives 



S(6, A, t) = ^G s ®il + ^G S ® S + S E ~ ^G S <E)J - ^G S ® S, (17) 
7 7 7 7 

which represents the SLE in the "gravitationally self-consistent" form that is numerically imple- 
mented in SELEN. Since the unknown S(w,t) also appears in the spatiotemporal convolutions 
at the right-hand side, the SLE is an integral e quation, which cannot be solved explicitly unless 
some drastic simplifying assumptions are made (jSpada and Stocchi l200fih . The SLE is a linear 
equation as long as shorelines are not allowed to migrate horizontally, i.e. if O (and consequently 
A ) are not time-dependent. Sea level variations are se nsitive to mantle rheology through G 



since this is determined by the viscoelastic LDCs h and k ( Spada . 2003; Spa da and Stocchil . l2006l ). 
Solutions of Eq. (|T7|) in special cases, discussed in detail by ISpada and Stocchil (|2006l ). are also 
available via SELEN . 
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3 How SELEN works 



Running SELEN requires a standard UNIX (including Linux and Mac OS X) environment and 
a Fortran 90 compiler. On Windows systems, SELEN can run within the Cygwin environment, 
freely downloadable from http://www.cygwin.com. A list of supported operating systems and For- 
tran compilers is given in the (simple text) configuration file config.dat. While the previous 
version of SELEN (1.0) was limited to the IBM XL commercial compiler, SELEN 2.9 compo- 
nents can be compiled both using the freely available g95 (http://www.g95.org) and gfortran 
(http://gcc.gnu.org/gfortran) compilers, or the commercial Intel Fortran compiler. Additional 
configurations for the operating system and compilers in SELEN 2.9 can be implemented by mod- 
ifying the setup progr am conf ig.f90. To ru n SELEN, the GMT (Generic Mapping Tools) public 
domain software of Wess el and Smith! (Il998h must be installed on the system. In SELEN 2.9 , 
the most computationally intensive portions of code have been parallelized with QpenMPl (2005) 
directives. The corresponding program units can therefore take advantage of multi-threading on 
modern CPUs, resulting in a substantial performance improvement. 

SELEN consists of several (independent) Fortran 90 program units, performing specific compu- 
tation steps. A list, with a short description of their purpose, is provided in Table [5j The user can T]5] 
specify run parameters in the configuration file config.dat; this file is parsed by the setup pro- 
gram config.f90, which checks parameters for consistency and creates a shell script (selen.sh) 
containing the compilation and execution commands for the SELEN components needed for the 
specific run. The configuration and execution process is transparent to the user and entirely han- 
dled by the makeselen. sh shell script. To launch a SELEN run, the user only needs to edit the 
config.dat file and execute the makeselen . sh script. This constitutes a significant improvement 
with respect to previous version SELEN 1.0, which lacked a flexible user interface. After each run, 
SELEN outputs are stored in subfolders of the output folder named depot -name, created in the 
working directory, where name is a 4-character label defined by the user in the config.dat file. 

In this way, outputs of different SELEN runs can be stored in the working directory. The 
structure of the output folder is described in Table [5j a copy of the config.dat file is also stored T]5] 
in the output folder for reference. 

A typical SELEN run can be subdivided in three main steps: i) pre-computing of all the needed 
functions, ii) numerical solution of the SLE and Hi) computation of geophysical and geodetic 
quantities requested by the user. 

During the first phase, SELEN performs the spherical harmonic (SH) expansions of the ocean 
function O and of the ice thickness variation /, and obtains the LDCs h, I and k for the selected 
Earth model (an Earth model is defined by its elas tic and rheo logical layering according to the 
conventions adopted for the TABOO program, see ISpadal 120031 ). For the numerical evaluation 
of surface integrals of Eq. @, SELEN defines an icosahedron-shaped pixelization of the sphere, 
displayed in Fig. [TJ The grid density is controlled by the resolution parameter R, defined in IQ] 
th e config.dat fil e; the number of grid pixels is given by N p = A0R(R, + 1) + 12. As discussed 
bv lTegmard lfl99fih. an accurate numerical integration is only ensured for £ max < \/3N p , where 
f-max is the maximum harmonic degree of the SH expansion. A discrete realization of the ocean 
function is practically obtained by selecting the grid pixels falling in "wet" and "dry" areas with 
the gmtselect utility from the GMT package (see Fig. [T|). To optimize computations involving 
SH functions, all the needed associated Legendre functions and trigonometric functions at grid 
pixels are pre-computed once and stored on disk for later re-use. This SH database can be re- 
used also across different runs, as long as parameters R and £ m ax are unchanged (values of R and 
(■max are defined in the config.dat file). Similarly, the SH expansion of the ice load, a relatively 
time-consuming task, can be re-used by setting the appropriate option in the configuration file 
config.dat. 

In the second phase, SELEN solves numerically the SLE. As discussed in Section [21 the SLE is 
an integral (implicit) equation, in which the unknown S is convolved in space and time with the 
sea level Green's functions. An iterative solution scheme is therefore needed, similar to the one 
employed in the solution of one-dimensional inhomogeneous integral Fredholm equation of the 
second kind. As a zero-order solution, SELEN approximates the unknown function S with the 
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eustatic sea level variation: = S E . By substituting in the right-hand side of Eq. (|TT)) . 
a new estimate is obtained. This process can be iterated any number of times, by using 
Eq. (IT71) to obtain the fc-th order solution from the (k — l)-th order solution S 1 ^ -1 - 1 , until 
a pre-determined convergence threshold is reached (i.e. until the ratio \S^ - |/|s« | 

becomes suffici ently small at all the grid pixels and time steps). According to the numerical tests 
carried out by ISpada and Stocchil (|2007l ) and to our experience, three iterations normally suffice 
for convergence. 

In the third and final phase, SELEN uses the solution of the SLE in the spectral domain 
to obtain predictions of geophysical and geodetic quantities by SH synthesis. According to the 
settings in the user configuration file config.dat, SELEN can compute present-day rates of 
deformation (in vertical and horizontal directions), geoid change and sea level change on global or 
regional scales, relative sealevel predictions at specific sites, sea level change rates at tide gauge 
stations, present-day velocities at user-specified sites, and present-day rates of variation of the 
Stokes coefficients of the Earth's gravity field. 

The execution time of SELEN scales with spatial and harmonic resolutions as t exe ~ N p Nh, 
where N p is the number of grid pixels and Nh — {£ ma x + ^){£max + 2)/2 is the number of SH 
functions with harmonic degree £ < £ max ', as a rough approximation, t exe ~ R 2 l 2 max ^ where R is 
the resolution parameter. The memory footprint of SELEN follows a similar scaling law. For high 
resolution runs (large values of £ ma x and consequently of R), disk I/O can become a considerable 
fraction of t exei and these relations may no longer be valid. On multi-core systems, SELEN can use 
multi-threading to reduce computation time (multi-threading was not implemented in previous 
version SELEN 1.0). To enable multi-threading, the corresponding option in config.dat must 
be set, specifying the number of threads that SELEN will create. The number of threads should 
be equal to the number of cores in the system; it can be smaller if the user wishes to leave some 
resources free for other tasks. 



4 A few case studies 



To illustrate the main feat ures of SELE N, we have configured the program for a GIA simulation 
based on model ICE-5G (jPeltien . 120041 ) . In the previous version (SELEN 1.0) this ice model, 
which is now widely employed in the GIA literature, was not accessible. The configuration file 
config.dat for this run (referred to as TEST run) and all the SELEN output files (text files and 
plots) are available in the folder depot -TEST that comes with the SELEN distribution package. 
The numerical results presented below are organized into two subsections. In the first (Section 
14. ip . the focus is on the Relative Sea Level (RSL) variations driven by GIA since the LGM, 21 kyrs 
ago. The second (Section I4.2[) is about the delayed viscoelastic effects that GIA is still producing 
today, and the focus is on the time variations of various geodetic quantities. Keeping separate the 
two time scales helps to simplify the presentation of the results. However, it should be kept in 
mind that the present-day geodetic variations depend on the whole history of deglaciation, and 
are thus extremely sensitive to mantle visco sity in spite o f the long time elapsed since the end of 
melting of major continental ice sheets (e.g. |Peltierll2004l ). 

The viscosity profile employed in the cases studies illustrated b elow is a thre e-layer volume- 
average of the original, multi-layered "VM2" profile introduced by iPeltierl f|2004h and associated 
with the ICE-5G model. Values of viscosity and other parameters that define the incompressible 
Earth model (hereafter referred as to VM2a) are summarized in Table [5] The isostatic relaxation 
spectrum determined by SELEN for VM2a is shown in Fig. [21 while F i g. El s hows the elastic and 
fluid LDCs. SELEN incorporates TABOO (|Spadal . 120031: ISpada et all \2004 as a subrou tine for 
computing these spectra, based on the Viscoelastic Normal Modes theory (jPeltierl . fl974l ). After 
execution, all the information about the LDCs (tables with ASCII data and plots, if requested) are 
made available to the user in folder depot-TEST/Love-Numbers-by-TABOO, including the viscous 
parts of the LDCs (not shown here). All the ensuing SELEN results are obtained using LDCs 
in the range of degrees 1 < l max < 128 and are expressed in the reference frame with origin the 
Earth center of mass (hence, differently than in SELEN 1.0, here we properly account for the 
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degree I = 1 LDCs). LDCs of harmonic degree I — vanish identically because of incompress- 
ibility, but they would not play any role even for an compressible Earth since the SLE includes 
explicitly the constraint of mass conservation (see Eq.[7]). The SLE is solved iteratively on the grid 
shown in Fig. [J and three iterations are performed to ensure convergence of the iteration scheme 



(jSpada and Stocchil . 120071 ) 



The spatial distribution of the ice thickness according to model ICE-5G at the LGM and at 
present time are shown in the SELEN plots of Fig. [4] The ICE-5G thickness data are obtained from E]4] 
the home page of Prof. W. R. Peltier (the version with 1° resolution is used). In folder ICE-MDDELS 
of the main SELEN directory, other ice models are available, including the previous ICE-X models 
developed by Prof. W. R. Peltier and co-workers (SELEN also incorporates individual components 
of the global ICE-5G model, which could be useful for regional studies). Using specific Fortran 
formats described in program conf ig.f90, the user can indeed introduce other ad hoc ice models 
according to specific purposes. In the specific case of ICE-5G, the spatiotemporal discretization 
has been slightly modified for the sake of co mputational convenience (details on the discretization 
scheme adopted in SELEN are available in lSpada and Stocchil 120071 ) . The original ICE-5G time 
grid has been converted into a uniformly spaced 1— kyr grid, and the elementary "rectangular" 
1° x 1° ice elements that define ICE-5G are converted, at run time, into equal-area and equal- 
thic kness "discs" that allow f or a straightforward SH decomposition because of their symmetry (see 
e.g. ISpada and Stocchill2006l) . Since in SELEN 2.9 fixed shorelines are assumed, the melting of the 



ice distributed across the continental shelf in Fig. @] (top) is not accompanied by a transgression 
of ocean water, which may imply a local error in sea level change predictions (this restriction will 
be relaxed in future releases of SELEN). The equivalent sea level (ESL), shown in Fig. [5l provides E(S] 
the time-history of the ice volume of ICE-5G since the LGM. After execution, all the ice model 
data are stored in depot-TEST/ICE5G. These include the spatial distribution of the ice masses at 
all times, the ESL function, the SH coefficients of the ice thickness and the SH reconstruction of 
the ice distribution (options for these computations are available in file config.dat). 

4.1 Past Relative Sea Level variations 

One of the purposes of SELEN is the modeling of sea level variations that occurred during the 
time period elapsed since the LGM, in consequence of the melting of late-Pleistocene ice sheets. 
There is a considerable amount of literature about this problem, and it is impossible to provide an 
exhaustive summary in this paper. A flavor of the enormous importance t hat past sea level varia- 



tions h ave on our unders tanding of present-day variations is given by e.g. lLambeck and Chappell 



(|200lf ) and iPeltierl (|2004l ) and references therein. 

In this section, we provide an illustration of the SELEN outputs regarding "Relative Sea Level" 
(RSL). This quantity is not directly obtained from the SLE (flTl) . which provides S(ui, t). According 
to the geological practice, RSL is defined as the difference between SL at a given epoch before 
present (t = tpp) and the present -time (t — t p ) value: 

RSL(w, t BP ) = SL(w, t BP ) - SL{u>, t p ), (18) 

which, using the definition of sea level given by Eq. ([1}, can be also written as 

RSL(w, t BP ) = 5(w, t BP ) - S{u, t p ), (19) 

showing that RSL can indeed be obtained from the solution of the SLE computed at two differ- 
ent times (we note that RSL does not depend on the choice of the remote time t r in Eq. @). 
All the following RSL computations have been obtained performing three iterations within the 
solution scheme of the SLE, which genera lly provides a sufficiently accurate solution (see the test 



computations in lSpada and Stocchil 2007) 



SELEN can be easily configured for producing RSL predictions at sites of interest to the user, 
for which observations of past sea levels can be obtained from the literature. Fig. [6j produced by F|6] 
SELEN, shows the locations of 392 sites for which radiocarbon-controlled RSL data are available, 
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according to the compilation of Tushingham and Peltier ( 19921 Il993). Information on the sites co 



ordinatcs and RSL observations, including their uncertainties, are contained in file sealevel.dat, 
obtained from the National Oceanic and Atmospheric Administration (NOOA) page [J The RSL 
observations in this file cover the last 15,000 years and are obtained from various sources, mainly 
based on geomorphologic and archaeological methods. In the scatterplot of Fig. all RSL obser- EJ7J 
vations from the input file sealevel . dat (top frame) are qualitatively compared with predictions 
obtained from SELEN using the settings in config.dat (bottom). The scatterplot gives a clear 
view of the temporal distribution of the RSL observations from this data collection, which are 
mainly relative to the last ~ 8, 000 years. The similarity of the two scatterplots in this figure 
clearly indicates that the ICE-5G model, in its SELEN implementation, broadly reproduces the 
RSL observations globally available. Of course, local misfits are possible, as we will discuss below 
for specific sites. A more rigorous global misfit analysis is also possible using SELEN, oriented 
to users interested to tackle an optimization problem in which the best fitting model parameters 
(i.e. the viscosity of mantle layers) are to be inverted from a specific set of RSL observations 
(the results of the misfit analysis for this run, not reproduced here, are available in the folder 
depot-TEST/rsl/rsl-misf it). 

Fig. [5] shows a RSL analysis, in which SELEN is programmed to produce RSL predictions at E[5] 
individual sites. Here we only show eight curves out of the 392 available in folder depot-TEST/rsl/ 
rsl-curves after the execution of SELEN, all relative to the case of the Hudson Bay. The 
process of GIA across the Hudson Bay has been the subject of various investigations in the past, 
because of the sensitivity shown by the uplift d ata on the rheological layering of the mantle 



( Mitrovica and Peltier . 1992 : Icianetti et al. . 2002 ). RSL observations from this region, shown by 
error bars in Fig. clearly indicate a monotonous sea level fall that is typical of sites belonging to 
previously ice-covered areas during the LGM. Model predictions reproduce the RSL observations 
very satisfactorily. This should not come as a surprise, since the two basic components of model 
ICE-5G (namely, the chronology of the late-Pleistocene ice sheets and the viscosity profile of 
the mantle) have been expressly designed to be st fit a global dataset of RSL observations, which 
includes these Hudson bay sites ( Peltier . 2004fl . The further example of Fig. [5] shows results for F(5] 



eight RSL sites in the far field of the previously ice-covered regions, where the history of sea 
level rise has been less influenced by the direct effect of ice melting. The shapes of the RSL curves 
in these regions clearly depart from those of the Hudson bay in Fig. [51 In some cases, they neatly 
indicate a monotonous sea level rise, but some remarkable exceptions are found, in which they 
show sea level high-stands or more complex features. Their trend is reasonably reproduced by 
our SELEN computations. 

In Fig. [10] we show the results of further runs, in which the program has been configured I fTOl 
to solve the SLE using various approximations. The ice model and the rheology are the same 
as in previous computations. This kind of analysis, which is performed here only for the sites 
of Richmond (Hudson bay) and of Merseyside (England) , is not meant to test the agreement of 
RSL predictions with observations. Rather, it can be useful to illustrate the role played by some 
of the basic physical ingredients of the SLE. The solid curves still show the "gravitationally self- 
consistent" (GSC) solution, in which the SLE is solved in the full form given by Eq. (fTT)) . In results 
shown by dotted RSL curves, only the elastic (EL) components of the LDCs are employed in the 
computations, hence neglecting the viscous components of the Maxwell rheology. The eustatic 
curve (EU, dashed), shows the RSL trend obtained assuming a rigid Earth and neglecting any 
gravitational interaction between the solid Earth, the ice masses, and the oceans. The eustatic 
RSL curve, which is simply expressed by Eq. [TU1 is the same for all RSL sites. Finally, the dash- 
dott ed curve (W) shows the solution of the SLE in the so-called "Woodward approximation" 
(see Spada and Stocchil l2007t ) . in which the Earth is assumed to be perfectly rigid (LDCs are 



h = I = k = 0), and only the gravitational attraction between the ice masses and the oceans is 
taken into account. From the results of Fig. [TTjl it is apparent that the RSL observations across 
the Hudson bay, where ice melting has produced a huge uplift, can only be explained invoking the 
"gravitationally self-consistent" solution. However, in the case of the far-field site Merseyside, the 



1 See http://wwwl. ncdc.noaa.gov/pub/data/paleo/paleocean/ relativeseajevel/. 
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GSC solution can be assimilated to any of the approximate solutions (including those pertaining 
to a rigid Earth, EU and W), at least for the last few kiloyears. 

The site-by-site analyses shown in Figs 191 and ITUl above are useful for a detailed reconstruction 
of the local RSL curves and direct comparison with observations. However, a global or regional 
visualization of the RSL variations can be more useful for a qualitative interpretation. For the 
purpose of visualization of the relative sea level variations associated with GIA, SELEN can be con- 
figured in two different ways. Both features were not directly accessible to the user in SELEN 1.0. 
The firs t, which is o r iented to global analyses, allows for the visualization of the so-called "Clark 



zones" (| Clark et all 119781 ). i.e. the regions of the globe in which the RSL curves shows similar 



patterns after the LGM. Numerical results and plots for this analysis, which are not reported here 
for space limitations, will be accessible for the TEST run in folder depot -TEST/rsl/rsl-zones/ 
after execution of SELEN. An example of the second possible RSL analysis, which is more oriented 
to regional investigations, is shown in Fig. [II] where RSL contour plots across the Mediterranean 
region for times 2, 6 and 10 kyrs BP are drawn. This area is the subj ect of considerable in- 



terest in view of the signifi cant amount of sea level indicators available ([Lambeck et all [2004; 



Lambcck and Purccll, 2005). Fig. [TT] provides a clear regional characterization of the process of 



sea level rise in the area. Although the RSL pattern shows a considerable complexity, it is appar- 
ent that the largest RSL excursions are predicted in the b ulk of the basin, in conseq uence of the 
broad subsidence caused by the effects of meltwater load (jStocchi and Spad"al 120071) . Data and 
maps for this kind of analysis are found in folder depot -TEST/rsl-contours/. 



En] 



4.2 Effects of GIA on present— day geodetic variations 

Three-dimensional movements of the solid Earth and gravity variations are still affected signif- 
icantly by the isostat i c dise quilibrium in response to the melting of late-Pleistocene ice sheets 
(see e.g. King et al. 2010fl . Therefore, it is sometimes necessary to evaluate quantitatively 



the amplitude of the GIA effects at present time, possibly in order to decontaminate geodetic 
quantities observed locally or globally. Since for loads of the size of major late-Pleistocene ice 
sheets the Maxwell relaxation time of the mantle is of the order of a few kilo-years (see e.g. 
Schubert et al. 2004), the rates of sea level change or other quantities associated with GIA can be 
effectively considered as constant through a decade to century time scale. Hereinafter, we will con- 
sider a few examples in which the SLE is solved for trends of GIA-related quantities. Since solvin 
the SLE implies a temporal discretization of all variables involved (see ISpada and Stocchil |200 
for details), in SELEN geodetic trends at a specific place ui are evaluated numerically as 



Q(uj,t p ) - Q(uj,t p - A) 
A 



(20) 



where A = 1 kyrs is the natural time step in SELEN, and Q here represents any of the fundamental 
geodetic quantities S (relative sea level variation), U (vertical disaplacement) and N (sea surface 
variation) that appear in Eq. [3] Of course, more accurate methods can be implemented by the 
user for the numerical evaluation of trends. Thought here we focus in the effects of melting of 
past ice sheets on present-day geodetic vari ations, with min o r modifications SELEN can be also 
employed for the study of the fingerprints ( Mitrovica et all 12001 ; Tamisiea et al. L 2001 ) o f the 
current terrestrial ice melting. Examples have been recently given by Sorensen et al.l ( 2011 ) and 
ISpada et al.l (2012) for the case study of the melting of the Greenland ice sheet. 

Present-day trends S, U and N, computed by means of Eq. (j2TJ)) . are shown in Figs [T2l Upl and 
ITU respectively, for our TEST run. Of course, these global maps are not independent of one another, 
since from the SLE (see Eq. [3]), we have e. g. N = S + U. To emphasize the regional variations 
of these GIA signals, the color tables span the range of ±1 mm/yr. The very complex pattern of 
sea level change S in Fig. [12] is the consequence of the delayed visco-elastic deformations of the 
solid surface of the Earth, of gravitational interactions between the solid Earth, the ice masses and 
the oceans, and of the complex shape of continents. For a rigid, non gravitating Earth we would 
observe 5 = ever ywhere on this map. The S pattern in Fig. [T2] has been discussed in a number 
of studies (see e.g. Mitrovica and Milne] 2002 ) . Some features are easily interpreted, such as the 
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broad areas of sea level fall across the formerly glaciated regions and the surrounding sea level rise 
corresponding to the collapsing forebulges. Some other features, such as th e equatorial region of 
sea le vel fall, have a more difficult interpretation in terms of ocean siphoning (jMitrovica and Milnel . 



2002) 



The vertical velocity (U) map in Fig. 1131 clearly shows a strong anti-correlation with S, and in 
particular it very neatly shows that, in the far field of previously glaciated areas, the continents 
are currently moving up relative to the Earth's center of mass at rates as large as 0.5 mm/yr. This 
is ultimately the consequence of the differential movements induced by meltwater loading acting 
on the oceans floor. Compared with S and U, the rate of sea surface variation N in Fig. Q3] is 
characterized by a smoother pattern globally, which manifests a relatively large energy content 
of low-degree harmonics. As seen from the Earth's center of mass, the sea surface is collapsing 
everywhere at rates between and 0.5 mm/yr, except in the formerly glaciated regions, where we 
observe an uplift that broadly follows the uplift of the solid surface of the Earth in Fig. Q21 but 
with a significantly smaller amplitude. 

In addition to global maps of trends S, U and N, SELEN can similarly produce regional 
analyses, which were not directly accessible in previous release. An example is shown in Fig. ll5l for 
the Mediterranean region. This figure shows very clearly how GIA is acting across this relatively 
small, mid-latitude basin: the region is subsiding (U < 0, middle frame), and though the sea 
surface is collapsing relative to the Earth's center of mass (N < 0, bottom), in the bulk of the 



Medit erranean relative sea level is rising (S > 0, top). As discussed in e. g. IStocchi and Spada 
(2009), meltwater loading, described by the second term on the right hand side of Eq. (fTTj) , is the 
major cause of subsidence. This is indicated by the shape of the contour lines of U, which are 
broadly following the coastlines. 

Tide gauge (TG) records hold a central role in all the assessments of the secular global mean sea 
level rise so far published (see Douglasfl997 , Spada and Galassl2012 and references therein). Since 



TGs measure the offset between sea surface and the solid Earth, they constitute the experimental 
realization of the SLE in its basic form, given by Eq. ([3]). Perhaps it is not fully appreciated by the 
climate cha nge community that the current estimates of a global sea l evel rise of 1.8 ± 0.1 mm/y r 
since 1880 ( Douglas! 1997 ). recently revised to 1.5 ± 0.1 mm/yr by ISoada and Galassl ( 20121 ). 



tightly depend on GIA modeling (hence, to a large extent, on the rheology of the solid Earth). 
In fact, GIA models consistent with RSL observations since the LGM are commonly employed 
to decontaminate TG trends obtained from long re cords, in orde r to fully highlight the effect of 



secular climate variations on sea level rise (see e. g. lPeltierll2001l) 



Within SELEN, the GIA component of sea level change at TGs is evaluated by 



r 



gta 
k 



S(u k ,t p ), (21) 



where t v is present time, the time derivative of sea level change is given by Eq. (|20p and lju = 
(9k, Xk) denote colatitude and longitude of the k—th TG. Since the Maxwell relaxa ti on tim e of 



the bulk of the mantle is of the order of a few kilo-years (see e.g. ISchubert et al.l 120041 ) . the 



rate (|2ip can be effectively considered as constant through the period of the instrumental TG 
record. In SELEN, a simple analysis can be performed in which predictions of r 9 k ia are computed 
at the TG locations. In the TEST run of SELEN, we consider the input data collected in file 
DATA/rlr-trends .txt, pertaining to the 1123 Revised Local Reference (RLR) TG stations listed 
by the Permanent Service for the Mean Sea Level (PSMSL) as of January 22, 2007 (SELEN maps 
for these data are shown in Fig. ll6[) . Of course, the user can input other data sets, not necessarily 
containing TG information. Sample outputs of this SELEN analysis are shown in Tables [5] and 
[5j The first shows results for TGs with records of at least 100 years of RLR observations, while 
the second shows a selection of Mediterranean TGs with at least 30 years of observations. The 
two tables are showing GIA predictions for the rate of sea level change r^ za , but also for the 
other fundamental geodetic quantities U and N (a new feature with respect to SELEN 1.0). All 
data and plots produced by SELEN for the TG analyses illustrated above are found in folder 
depot-TEST/tgauges. 

As a final example of a possible geodetic application available within SELEN, we consider the 
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time-variations of the gravity potential harmonic coefficients. This analysis can be useful, for ex- 
ample, to evaluate the effects of GIA on the gravity variations observed by the NASA/DLR Gravit y 
Recovery and Climate Experiment (GRACE) satellites (see the overview of lTaplev et ajj I l2009h . 



According to SELEN conventions, the GIA-induced total variation of the gravity field potential, 
evaluated at the Earth's surface is: 

p p 

S(a,w) = -f ]T E QzrrMu), (22) 

t=2 m=-t 

where T is Newton's constant, M e and a are the mass and the average radius of the Earth, £ max 
is the maximum harmonic degree of the analysis (for run TEST, £ max = 128), <fr^ TO are the gravity 
coefficients in complex form (these obey <fr^_ TO =<I?^ TO ) and D^m(^) are the complex 47r-normalized 
SH functions of harmonic degree £ and order m: 



yiM = J(2£+l) ( f. , m) *: P em (cosO)e m \ (23) 
y (t + my. 

where Pp m (cos6) are the associated Legendre functions. Note that in Eq. (|2"2"j) . terms of harmonic 
degree £ = are not included since the total mass of the Earth is conserved, and those of degree 
1=1 vanish since the origin of the reference frames is assumed to coincide with the Earth's center 
of mass. It is sometimes more convenient to transform Eq. (|22j) into an equivalent expansion 
involving real coefficients. For instance, according to the conventions usually adopted in gravimetry 
(including GRACE), we can equivalently write 

$(a, w) = V* (ctm cosmA + sp m sinmA) P| m (cos6>), (24) 

t=2 m=0 

where Pg m (cos6) denotes the "fully normalized" associated Legendre functions, not including 
the Condon-Shortley phase (—1)™. After some straightforward algebra, for the variations of 
the Stokes coefficients one obtains (c£ m ,sp m ) = (— l) m \/2 — So, n $| m - Program stokes. f 90 
provides, degree-by-degree, the time derivatives cp m and st m , which can be directly compared 
with trends obtained from monthly GRACE observations (derivatives are numerically evaluated 
following Eq. |2"0)) in order to assess the GIA effects. Fig. [T7] shows the result obtained for run 
TEST, for harm onics up to degree £ ma x = 9. A similar analysis has been performed recently by 
ISorensenl (|2010l h using various global ice models from the literature (and available in the ice models 
directory ICE-MODELS of SELEN). 



5 Conclusions 

We have described the basic features and the essential theory background of SELEN, an open 
source Fortran 90 program that solves numerically the SLE. Version 2.9, which constitutes a sig- 
nificant improvement of previous version SELEN 1.0, has been developed to respond to the requests 
of colleagues working on various aspects of GIA, who did not have access to a SLE solver, and for 
te aching purposes. SELEN , which is the numerical implementation of the SLE theory presented 
bv lFarrelland Clarki ^97^), can be used to efficiently simulate the spatiotemporal variations of 



several geophysical and geodetic quantities involved in the GIA process. These include relaxation 
spectra and LDCs computed using user-supplied rheological profiles, the spatial distribution of 
the continental ice sheets according to several built-in GIA models, predictions of relative sea 
level variations at specific sites from which observations are available since the LGM, global and 
regional maps of present-day rates of variations of GIA-induced sea level change and displace- 
ments, local predictions for geodetic variations (including vertical and horizontal movements), 
and time-variations of harmonic coefficients of the gravity field in response to GIA. Refined im- 
plementations of SELEN 1.0 have been recently employed in various geodetic contexts, ranging 
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from the interpretation of GPS data on a regional scale in Greenland ([Nielsen et all l2012i ) . to the 



correction of GRACE observations aimed at estim ating the mass variation in the Mediterranean 
and in the Black Sea (IFenoglio-Marc et all l2012f) . and to the study of the effects of GIA in West 
Antarctica (|Groh et all l2012h and Greenland (|Ewert et all 120121 ISpada et ail . 120121) . 

SELEN is written in standard Fortran 90 and takes advantage of an OpenMP multi-threaded 
parallelism. It can be installed and compiled on Mac OS, Linux and Windows platforms. Starting 
from the examples illustrated in this manuscript, the users can customize SELEN in order to 
solve ad-hoc GIA problems or to include new functionalities within the code. According to our 
experience with the students of the 2nd Training School on Glacial Isostatic Adjustment (GIA) 
Modeling, where SELEN has been employed for a tutorial introduction to the physics and the 
phenomenology of GIA, the program may constitute a valuable instrument for scientists at their 
first approach to this problem. New features such the simulation of the horizontal migration 
of shorelines in response to sea level change and the realization of the feedback between sea 
level variations and rotational fluctuations, are under way and will be available with the upcoming 
version of SELEN (v. 3.0). 
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config.dat SELEN configuration file 

makeselen.sh Bash script that executes SELEN 



OT^T I'M T , 1 j 

SELEN Input data: 




wdir/DATA 


Various input data for RSL analyses and more 


wdir /ICE-MODELS 


Ice models data 


wdir/VSC 


Collection of viscosity profiles 


Fortran 90 units: 




conf ig. f 90 


Interpreter ol the conhguration hie 


esl . f 90 


Equivalent Sea Level 


geo.f90 


Time variations of geodetic quantities 


gmaps . f 90 


Synthesis ol geodetic quantities on global maps 


harmonics . f 90 


Include hie with Various SH tools and utilities 


ms.f90 


GMI multi-segment hies horn ice data 


of _dv.f 90 


Degree variance ol the Ocean function 


px . f 90 


Pixehzation tools (including the legmark algorithm) 


px_rec . f 90 


Reorganizes the pixelization data 


rec_ice . tyu 


SH reconstruction of the ice thickness 


rec_of .f 90 


SH reconstruction of the Ocean Function 


rmaps . f 90 


Synthesis of geodetic quantities on regional maps 


rsl.f90 


Relative Sea Level curves 


rsl_zones . f 90 


Geometry of the Relative Sea Level "Clark's zones" 


rslc.f90 


Relative Sealevel Contour lines for regional analyses 


sh.f90 


SHs at the grid pixels 


sh_of .f90 


SH coefficients for the Ocean Function 


shjrsl.f90 


SHs at the RSL sites 


shjrslc.f90 


SHs for regional analysis and RSL contours 


sh_tgauges . f 90 


SHs at the TG sites 


shape_f actors . f 90 


"Shape factors" for the ice elements 


shice.f90 


SH decomposition of the ice model 


shtools.f90 


An SHTOOLS interface for the SH analysis 


sle.f90 


The SLE solver 


stokes. f 90 


Variations of the Stokes coefficients of the gravity field 


tgauges . f 90 


Present-day rate of sea level change at the TG sites 


tb.F90 


The TABOO code 


wnw.f90 


Numerical test for the SH orthogonality 



Table 1: Some of the files contained in the SELEN package. Here wdir indicates the working 
directory of SELEN. 
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Folder : 



content : 



ICE5G 

ICE5G/esl 
ICE5G/original 
ICE5G/reconstructed 
ICE5G/sh 

Love-Numbers-by-TABDO 

TABOO 

bin 

geod 

geod/3dmaps 
geod/sites 

gmaps 

log 

of 

of /degree_variance 

px 

rmaps 

rmaps/Antarctica 



Data about the ice model 
ESL 

Ice thickness data 

SH reconstruction of ice thickness data 
SH coefficients of the ice model 
Love numbers data 
TABOO input hies 
An empty folder 
Predictions at geodetic sites 
Regional maps (in progress) 
Geodetic predictions at specific sites 
Global geodetic rates (data and plots) 
Log files of SELEN and TABOO 
(with summary of Earth parameters) 
Ocean function data and plots 
Ocean function degree variance 
Various pixelization data and plots 
Regional geodetic rates (data and maps) 
rates for Antarctica 



rmaps/North_America 

rsl 

rsl-contours 

rsl-curves 

rsl-misf it 

rsl-scplot 

rsl-sites 

rsl-table 

rsl-zones 
stokes 
tgauges 

tgauges-predictions 

tgauges-scplots 

tgauges-sites 

wnw 



rates for North America 

Relative Sea Level (RSL) data folder 

RSL contour plot 

RSL curves at specific sites 

Misfit between RSL data and predictions 

Scatterplot of RSL data 

Data and plots regarding RSL sites 

Summary table of RSL data and predictions 

RSL zones 

Stokes coefficients data 
Tide gauges (TGs) 
Predictions at TGs 
TG data scatterplot 
Maps of TG sites 

"Window test" for the ocean function 



Table 2: Organization of the SELEN outputs in directory depot -TEST. SELEN outputs consist of 
various plain text, postscript, and pdf files. 
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RES= 44, N= 75692 



BS3 201 1 Dec 20 08:55^34] SELEN 2.9 



Figure 1: Pixelization of the sphere following the icosahedron-based method proposed bv lTegmark 



(199<§) for astrophysical applications, implemented in SELEN. This geometrical tool provides a 
natural set of Gauss points on the sphere and allows for a straightforward computation of integrals 
involving SH functions. The number of pixels in the grid is N p — 40R(R + 1) + 12, where R is the 
resolution parameter. Here R = 44 (the corresponding grid spacin g is ~ 45 km). The constraint 
N p > £max/3, which ensures an optimal integration on the sphere ( TegmarkL Il996 ) is largely met 



(in the TEST run, £ max = 128). Data about the grid, including spherical coordinates of "wet" 
(oceanic, blue) and "dry" (continental, green) pixels and postscript figures are found in folder 
depot-TEST/ px. In SELEN, wet pi xels are separated by dry pixels using the GMT program 
gmtselect (jWessel and Smith! Il998l) . By default, SELEN employs the full resolution coastlines 
of GMT (-Df), and dry (wet) pixels are selected using option -Ns/k/s/k/s (-Nk/s/k/s/k) of 
gmtselect. 
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Layer 


Radius 


Density 


Shear modulus 


Viscosity 


Gravity 




(km) 


(kg m" 3 ) 


(xlO 11 Pa) 


(xl0 21 Pa s) 


(m s- 2 ) 


Lithosphere 


6281-6371 


4120 


0.73 


oo 


9.707 


Shallow upper mantle 


5951-6281 


4120 


0.95 


0.5 


9.672 


Transition zone 


5701-5951 


4220 


1.10 


0.5 


9.571 


Lower mantle 


3480-5701 


4508 


2.00 


2.7 


9.505 


Core 


0-3480 


10925 


0.00 





10.622 



Table 3: Model parameters for model VM2a, employed for the TEST run of SELEN. The whole 
library of models available within SELEN is accessible in the Fortran unit tb.F90. 
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Earth model: NV=3 C0DE=2 
Viscosity profile: /2. 7 0.5 0.5/*1E21 Pa.s 
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m EHflM 201 1 Oct 1 8 1 6:35^43] SELEN 2.9 
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Figure 2: Isostatic relaxation spectrum for model VM2a (see Table O, showing the relaxation 
times as a function of harmo nic degree I = n in the range 1 < I < 128. The physical meaning of 
the spectrum is discussed by ISpada (|2003l) . 
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Earth model: NV=3 C0DE=2 
Viscosity profile: /2. 7 0.5 0.5/*1E21 Pa.s 
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Figure 3: Elastic (left) and fluid (right) values of the LDCs h (associated to vertical displace- 
ment), I (horizontal displacement) and k (incremental gravitational potential) for model VM2a 
(see Table [5]), as a function of harmonic deg ree £. Note that h is normalized by (21 +1). For the 
definition of the LDCs see e.g. ISpadal ( 20031 ). 
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Figure 4: Ice thickness T(lu, t) at the LGM (21 kyrs ago, top) and at during the most recent time in - 
crement, between 1 kyr BP and present time (bottom), according to model ICE-5G (jPeltieii 120041 ) . 
The time-history of the Equivalent Sea Level for this model is shown in Fig. [5j Maps for all the 
other time steps between LGM and present are available in folder depot-TEST/ICE5G/original. 
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Figure 5: Equivalent Sea Level for model ICE-5G ([Peltier . At a given time t before present, 

ESL(t) = (pi/ p w ){Vi{t) — Vi(t p ))/A ), where Vi{t) is the ice volume, Vi(t p ) is present day volume, 
and A Q is the area of the oceans surface. Hence, according to Eq. ([TO]), the plot of ESL mirrors 
that of S E . The total ESL variation (~ 127 m) represents the difference between ESL at the LGM 
and the present day value. 
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Figure 6: Geographical distribution of the 392 sites in file sealevel.dat, from which information 
about the history of RSL during the last <~ 15,000 years is available. This plot and more data 
about these sites are available from depot-TEST/rsl/rsl-sites. 
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RSL data from file: ./DATA/sealevel.dat 




61^2011 Nov 29 13:18:17 SELEN 3.2 



Figure 7: Scatterplot showing RS L observations (top) from sites of the compilation of 



Tushingham and Peltier! (|1992l I1993I ). RSL predictions, obtained solving the SLE in our TEST 



run of SELEN, are shown in the bottom frame. 
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/01 57.0 -77.0 1 RICHMOND Gl 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: 12.7 0.5 0.5/ 
-LMAX =128-RES .44 -ITER . 
-NV =3 -CODE =2 




102 53.0 -79.0 7 JAMES BAY QL- 
Ice model: ICE5G 
Repository label: FIG1 



-NV =3 -CODE =2 
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103 55.0 -82.5 7 C. HENRIETTA 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity prolile: /2.7 0.5 0.5/ 
-LMAX =128-RES .44 -ITER .3 - 
-NV=3 -CODE =2 
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104 58.0 -94.0 7 CHURCHILL W 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: /2.7 0.5 0.5/ 
-LMAX .128-RES =44 -ITER .3 - 
-NV =3 -CODE =2 
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105 64.5 -95.0 4 KEEWATIN NW 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: /2.7 0.5 0.5/ 
-LMAX .128-RES .44 -ITER . 
-NV =3 -CODE =2 
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106 59.8 -80.3 10 OTTAWA IS. N 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: /2.7 0.5 0.5/ 
-LMAX .128-RES =44 -ITER .3 - 
-NV =3 -CODE =2 
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107 64.5 -84.0 7 SOUTHAMPTO 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: /2.7 0.5 0.5/ 
-LMAX =128-RES .44 -ITER .3 - 
-NV=3 -CODE =2 
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108 62.0 -75.0 7 UNGAVA PEN. 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: /2.7 0.5 0.5/ 
-LMAX .1 28-RES .44 -ITER .3 - 
-NV =3 -CODE =2 
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Figure 8: RSL observations (with error bars) pertaining to the eight sites of Hudson bay 
in file sealevel.dat, compared with SELEN predictions (solid curves). Basic parameters 
for the TEST run are summarized in each frame. Postscript and PDF figures are located 
depot-TEST/rsl/rsl-curves/ps and depot-TEST/rsl/rsl-curves/pdf , respectively. 
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5571.6 103.4 7 SINGAPORE 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: 12.7 0.5 0.5/ 
-LMAX =128-RES .44 -ITER . 
NV =3 -CODE =2 



418 42.0 12.0 1 CIVITAVECCHIf_ 
Ice model: ICE5G 
Repository label: FIG1 
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528 -3.0 10.0 4 GABON 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: /2.7 0.5 0.5/ 
-LMAX =128-RES .44 -ITER .3 
-NV=3 -CODE =2 





358 32.0 -65.0 10 BERMUDA 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: 12.1 0.5 0.5/ 
-LMAX .128-RES .44 -ITER ..' 
-NV =3 -CODE =2 
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255 53.5 -3.2 8 MERSEYSIDE E 
Ice model: ICE5G 
Repository label: FIG1 
Viscosity profile: /2.7 0.5 0.5/ 
-LMAX .128-RES .44 -ITER .3 - 
-NV=3 -CODE =2 
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Figure 9: RSL curves at eight miscellanea sites located in the far field of the former ice sheets: 
Singapore, Civitavecchia (Italy), Gabon, NW Gulf of Mexico, Bermuda, Perth W. Australia, 
Merseyside (England), and Palmer peninsula (Antarctica). All the RSL predictions for the TEST 
run are found in depot-TEST/rsl/rsl-curves. 
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Figure 11: Example of RSL contour plot for the Mediterranean region, at times 2, 6, and 10 kyrs 
BP. These correspond to three different configurations of file DATA/rsl-region.dat, all based on 
the SELEN parameters employed in the TEST run (ice model ICE-5G and rheological parameters 
in Table EJ. 
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Rate of sea level change today 

-Ice model: ICE5G -Viscosity profile: /2.7 0.5 0.5/ 




-LMAX=128 -RES=44 -NV=3 -CODE=2 -MODE=1 -ITER=3 



Rate of vertical displacement today 

-Ice model: ICE5G -Viscosity profile: /2.7 0.5 0.5/ 




-LMAX=128 -RES=44 -NV=3 -CODE=2 -MODE=1 -ITER=3 




-LMAX=128 -RES=44 -NV=3 -C0DE=2 -M0DE=1 -ITER=3 




Figure 15: Regional analysis showing rates of sea level change (S, top), of vertical uplift (U, 
middle) and sea surface variation (AT, bottom) across the Mediterranean region, for run TEST. 
Data and plots for this analysis are found in folder depot-TEST/rmaps. 
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all tide gauge stations more than 30 years of data 

N= 1123 N= 460 




0" 180" 0" 0' 180" 0" 



EMBM2011 Nov 29 10:40:56 SELEN 2.9 

Figure 16: Geographical distribution of the TGs considered in run TEST, according to the number of 
RLR data available from each station. Data and plots for this analysis are available in depot-TEST/ 
tgauges/ tgauges-sites after the execution of SELEN. 
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Station Name valid yearly time range trend error r 9 ™ N U 





records 


year- 


-year 




mm 


/yr 




mm 


/yr 


mm 


/yr 


mm 


/yr 


Kungholmsfort 


118 


1887- 


-2005 


-0. 


.04 


0. 


.13 


-1. 


.48 


+0. 


.03 


+1. 


.52 


Olands Norra Uddc 


119 


1887- 


-2005 


-1 


.14 


0. 


.14 


-2. 


.60 


+0. 


13 


+2. 


.76 


Landsort 


119 


1887- 


-2005 


-2. 


.91 


0, 


.14 


-4. 


.18 


+0. 


.23 


+4. 


.44 


Nedre Sodertaje 


102 


1869- 


-1970 


-3 


.44 


0, 


.19 


-4. 


.73 


+0. 


.27 


+5. 


.02 


Stockholm 


117 


1889- 


-2005 


-3 


.90 


0. 


.15 


-4. 


.89 


+0. 


.28 


+5. 


.20 


Rat an 


112 


1892- 


-2005 


-7. 


.81 


0. 


.18 


-8. 


.95 


+0. 


.53 


+4. 


.49 


Helsinki 


124 


1879- 


-2004 


-2 


.43 


0. 


.16 


-4. 


.23 


+0. 


.29 


+4. 


.55 


Swinoujscie 


180 


1811- 


-1999 


+0. 


.82 


0, 


.06 


-0. 


.43 


-0. 


.08 


+0. 


.35 


Warnemunde 2 


149 


1856- 


-2005 


+1. 


.20 


0. 


.07 


-0. 


.36 


-0. 


.09 


+0. 


.27 


Wismar 2 


156 


1849 


-2005 


+1. 


.39 


0. 


.06 


-0. 


.26 


-0. 


.11 


+0. 


.15 


Kobenhavn 


104 


1889- 


-2002 


+0 


.39 


0. 


.11 


-0. 


.99 


-0. 


.02 


+0. 


.98 


Fredericia 


105 


1890- 


-2002 


+1 


.01 


0. 


.08 


-0. 


.57 


-0. 


.06 


+0. 


.52 


Aarhus 


101 


1889- 


-2002 


+0 


.55 


0. 


.08 


-0. 


.95 


-0. 


.03 


+0. 


.93 


Esbjerg 


103 


1890- 


-2002 


+1 


.18 


0. 


.14 


-0. 


.35 


-0. 


.08 


+0. 


.27 


Cuxhaven 2 


160 


1843- 


-2002 


+2 


.44 





.09 


-0. 


.01 


-0. 


.14 


-0. 


.12 


Aberdeen II 


103 


1862- 


-1965 


+0 


.58 


0. 


.10 


-0. 


.43 


-0. 


.17 


+0. 


.28 


Marseille 


108 


1886- 


-2004 


+1 


.28 


0. 


.08 


+0. 


.08 


-0. 


30 


-0. 


.39 


Poti 


121 


1874- 


-2004 


+6 


.55 


0. 


.15 


+0. 


03 


-0. 


.28 


-0. 


.31 


Batumi 


106 


1882- 


-2005 


+1. 


.78 


0. 


.20 


+0. 


.00 


-0. 


.28 


-0. 


.29 


Mumbay/Bombay 


111 


1878- 


-1993 


+0 


.78 


0. 


.09 


-0. 


.17 


-0. 


.27 


-0. 


.10 


Sydney, F. Denison 


108 


1886- 


-1993 


+0. 


.59 


0. 


.09 


-0. 


.32 


-0. 


.16 


+0. 


.15 


Honolulu 


101 


1905- 


-2005 


+1 


.48 


0. 


.12 


-0. 


.16 


-0. 


.35 


-0. 


.18 


Seattle 


106 


1899- 


-2005 


+2. 


.06 


0. 


.10 


+0. 


.71 


-0. 


.02 


-0. 


.73 


San Francisco 


151 


1855- 


-2005 


+1 


.43 


0. 


.08 


+0. 


.68 


-0. 


32 


-0. 


.99 


Baltimore 


102 


1903- 


-2005 


+3 


.11 


0. 


.10 


+1. 


.15 


-0. 


.10 


-1. 


.25 


New York 


131 


1856- 


-2005 


+2 


.76 


0. 


.06 


+1. 


13 


-0. 


.05 


-1. 


.19 



Table 4: Present-day trends of GIA-induced sea level change, sea surface variation, and vertical 
velocity at RLR PSMSL TGs with more than 100 valid years in the record series. With and 
<7fc we denote the best-fit secular rates and their uncertainties. This is an excerpt of the ASCII 
table depot-TEST/ tauges/ tgauges-predictions/ ptidegauges. 
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Station Name 


valid yearly 


time 


range 


trend 


error Ok 


r k 


N 


U 




records 


year- 


-year 


mm/yr 


mm/yr 


mm/yr 


mm/yr 


Algeciras 


34 


1944- 


-2001 


+0.46 


0.25 


-0.21 


-0.33 


-0.11 


Tarifa 


46 


1944- 


-2001 


+0.02 


0.39 


-0.21 


-0.33 


-0.12 


Malaga 


41 


1994- 


-2001 


+2.41 


0.43 


-0.23 


-0.32 


-0.09 


Alicante II 


30 


1960- 


-1995 


-0.85 


0.27 


-0.07 


-0.32 


-0.25 


Marseille 


108 


1886- 


-2004 


+1.28 


0.08 


+0.08 


-0.30 


-0.39 


Gcnova 


78 


1884- 


-1992 


+1.20 


0.07 


+0.07 


-0.29 


-0.36 


Venezia (S. Stefano) 


45 


1872- 


-1919 


+2.55 


0.42 


+0.04 


-0.27 


-0.31 


Venezia (P. Salute) 


82 


1909- 


-2000 


+2.39 


0.16 


+0.04 


-0.27 


-0.31 


Trieste 


96 


1905- 


-2006 


+1.17 


0.12 


+0.04 


-0.27 


-0.31 


Rovinj 


48 


1956- 


-2004 


+0.53 


0.29 


+0.06 


-0.27 


-0.34 


Bakar 


62 


1930- 


-2004 


+0.97 


0.23 


+0.05 


-0.27 


-0.32 


Split Rt Marjana 


50 


1953- 


-2004 


+0.60 


0.30 


+0.10 


-0.30 


-0.39 


Split Harbour 


50 


1955- 


-2004 


+0.33 


0.30 


+0.10 


-0.29 


-0.39 


Ceuta 


51 


1945- 


-2005 


+0.38 


0.23 


-0.21 


-0.33 


-0.11 



Table 5: Present-day trends of GIA-induced sea level change, sea surface variation, and vertical 
velocity at RLR Mediterranean PSMSL TGs with more than 30 valid years in the record series. 
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-Ice model: ICE5G -Viscosity profile: 12.7 0.5 0.5/ -LMAX=128 -RES=44 -NV=3 -CODE=2 -MODE=1 -ITER=3 
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1.0 



0.5 

x 
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2011 Nov 29 12:44:33 SELEN 3.2 
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j(l,m)=l(l+1)/2+m+1 



Figure 17: Time-derivatives of the Stokes coefficients of the Earth's gravity field associated with 
GIA, as a function of the generalized harmonic degree j = £(£ + l)/2 + m + 1 for 2 < t < 9 
and < m < i. Output data for this analysis are found in folder depot-TEST/stokes after the 
execution of SELEN. 



34 



References 



Cianetti S, Giunchi C, Spada G (2002) Mantle viscosity beneath the Hudson Bay: an inversion 
based on the Metropolis algorithm. J. Gcophys. Res. 107(B12):2352. doi:10.1029/2001JB000585 

Clark JA, Farrell WE, Peltier WR (1978) Global changes in postglacial sea level: a numerical 
calculation. Quat. Res. 9:265-287. 

Douglas BC (1997) Global sea-level rise: a redetermination. Surv. Geophys. 18:279-292. 

Ewert H, Groh A, Dietrich R (2012) Volume and mass changes of the Greenland ice 
sheet inferred from ICESat and GRACE. Journal of Geodynamics 59-60:111-123. doi: 
10.1016/j.jog.2011.06.003 

Farrell WE, Clark JA (1976) On postglacial sea-level. Geophys. J. R. Astron. Soc. 46:647-667. 

Fenoglio-Marc L, Rietbroek R, Grayek S, Becker M, Kusche J, Stanev E (2012) Water mass 
variation in the Mediterranean and Black Seas. Journal of Geodynamics 59-60:168-182. 

Groh A, Ewert H, Scheinert M, Fritsche M, Riilkel A, Richtcr A, Rosenau R, Dietrich R (2012) An 
investigation of Glacial Isostatic Adjustment over the Amundsen Sea sector, West Antarctica 
Global and Planetary Change 98-99:45-53. 10.1016/j.gloplacha.2012.08.001 

King MA, Altamimi Z, Boehm J, Bos M, Dach R, Elosegui P, Fund F, Hcrnandez-Pajares M, 
Lavallee D, Mendes Cerveira PJ, Penna N, Riva REM, Steigenberger P, Van Dam T, Vittuari 
L, Williams S, Willis P (2012) Improved constraints on models of Glacial Isostatic Adjustment: 
a review of the contribution of ground-based geodetic observations. Surv. Geophys. 31:465-507. 
doi:10.1007/sl0712-010-9100-4 

Lambeck K, Chappell J (2001) Sea Level Change Through the Last Glacial Cycle. Science 292:679- 
686. 

Lambeck K, Antonioli F, Purcell A, Silenzi S (2004) Sea-level change along the Italian coast from 
the past 10,000 years. Quat. Sci. Rev. 23:1567-1598. 

Lambeck K, Purcell A (2005) Sea-level change in the Mediterranean Sea since the LGM: model 
predictions for tectonically stable areas. Quat. Sci. Rev. 24:1969-1988. 

Milne GA, Mitrovica JX (1998) Postglacial sea-level change on a rotating Earth. Geophys. J. 
Int. 133:1-10. 

Mitrovica JX, Peltier WR (1991) On post-glacial geoid subsidence over the equatorial ocean. J. 
Geophys. Res. 96:20,053-20,071. 

Mitrovica JX, Davis JL, Shapiro II (1994) A spectral formalism for computing three-dimensional 
deformations due to surface loads. J. Geophys. Res. 99:7057-7073. 

Mitrovica JX, Peltier WR (1992) Constraints on mantle viscosity from relative sealevel variations 
in Hudson Bay. Gcophys. Res. Lett. 19:1185-1188. 

Mitrovica JX, Tamisiea ME, Davis JL, Milne GA (2001) Recent mass balance of polar ice sheets 
inferred from patterns of global sea-level change, Nature 409:1026-1029. 

Mitrovica JX, Milne GA (2002) On the origin of Late-Holocene highstands within equatorial ocean 
basins. Quat. Sci. Rev. 21:2179-2190. 

Nielsen K, Khan SA, Korsgaardb NJ, KH Kjasrb, Wahr J, Bcvis M, Stearns LA, Timm LH (2012) 
Crustal uplift due to ice mass variability on Upernavik Isstr0m, west Greenland Earth and 
Planetary Science Letters 353-354:182-189. doi:10.1016/j.epsl.2012.08.024 



35 



OpenMP (2005) OpenMP Application Program Interface, Version 2.5. OpenMP Architecture Re- 
view Board, [http://www.openmp.org/mpdocuments/spec25.pdf (last accessed 2011). 



Peltier WR (1974) The impulse response of a Maxwell earth. Rev. Geophys. Space Phys. 12:649- 
669. 

Peltier WR (2001) Global glacial isostatic adjustment and modern instrumental records of relative 
sea level history. In: Sea Level Rise: History and Consequences. Academic Press, San Diego, 
pp. 65-95. 

Peltier WR (2004) Global glacial isostasy and the surface of the Ice- Age Earth: the ICE-5G(VM2) 
model and GRACE. Annu. Rev. Earth PI. Sc. 32:111-149. 

Schubert G, Turcotte D, Olson, P (2004) Mantle convection in the Earth and Planets. Cambridge 
University Press, Cambridge. 

Solomon S, Qin D, Manning M, Chen Z, Marquis M, Averyt KB, Tignor M, Miller HL (2007) 
IPCC Climate change 2007: the physical science basis. Cambridge University Press, Cambridge. 

S0rensen LS (2010) Changes of the Greenland ice sheet derived from ICESat and GRACE data. 
PhD Thesis, University of Copenhagen. 

S0rensen LS, Simonsen SB, Nielsen K, Lucas-Picher P, Spada G, Adalgeirsdottir G, Forsberg R, 
Hvidberg, CS (2011) Mass balance of the Greenland ice sheet - a study of ICESat data, surface 
density and firn compaction modelling, The Cryosphere 5:173-186. doi:10.5194/tc-5-173-2011 

Spada G (2003) The theory behind TABOO. Samizdat Press, Golden, Colorado. 

Spada G, Antonioli A, Boschi L, Cianetti S, Galvani G, Giunchi C, Perniola B, Piana Agostinetti 
N, Piersanti A, Stocchi P (2004) Modeling Earth's post-glacial rebound. Eos. Trans. AGU 85:62- 
64. 

Spada G, Stocchi P (2006) The Sea Level Equation, Theory and Numerical Examples. Aracne, 
Roma. 

Spada G, Stocchi P (2007) SELEN: a Fortran 90 program for solving the "Sea Level Equation", 
Comput. and Geosci. 33(4):538-562. doi: 10.1016/j.cageo.2006.08.006 

Spada G, Barletta VR, Klemann V, Riva REM, Martinec Z, Gasperini P, Lund B, Wolf D, Ver- 
meersen LLA, King M (2011) A benchmark study for glacial-isostatic adjustment codes. Geo- 
phys. J. Int. 185:106-132. doi:10.1111/j.l365-246X.2011.04952.x 

Spada G, Ruggieri G, Sorensen LS, Nielsen K, Melini D., Colleoni F (2012) Greenland uplift 
and regional sea level changes from ICESat observations and GIA modelling, Geophys. J. 
Int. 189:1457-1474. doi: 10.1111/j.l365-246X.2012.05443.x 

Spada G, Galassi G (2012) New estimates of secular sea level rise from tide gauge data and GIA 
modelling. Geophys. J. Int. in press. doi:10.1111/j.l365-246X.2012.05663.x 

Stocchi P, Spada G (2007) Glacio and hydro-isostasy in the Mediterranean Sea: Clark's zones 
and role of remote ice sheets. Ann. Geophys. 50(6):741-761. 

Stocchi P, Spada G (2009) Influence of glacial isostatic adjustment upon current sea level variations 
in the Mediterranean. Tectonophysics 474:56-68. doi:10.1016/j.tccto.2009. 01.003 

Tamisiea ME, Mitrovica JX, Milne GA, Davis JL (2001) Global geoid and sea level 
changes due to present-day ice mass fluctuations, J. Geophys. Res. 106(B12):30, 849-30, 863. 
doi:10.1029/2000JB000011 



36 



Taplcy B, Bcttadpur S, Watkins M, Reigber C (2004) The gravity recovery and climate experiment: 
Mission overview and early results. Geophys. Res. Lett. 31(9):L09607. 

Tegmark M (1996) An icosahedron-based method for pixelizing the celestial sphere. ApJ Letters 
470:L81-L84. 

Tushingham AM, Peltier WR (1991) ICE-3G - A New Global Model of Late Pleistocene Deglacia- 
tion Based Upon Geophysical Predictions of Post-Glacial Relative Sea Level Change. J. Geo- 
phys. Res. 96:4497-4523. 

Tushingham AM, Peltier WR (1992) Validation of the ICE-3G model of Wiirm-Winsconsin 
deglaciation using a global data base of relative sealevel histories. J. Geophys. Res. 97:3285-3304. 

Tushingham AM, Peltier WR (1993) Relative Sea Level Database. IGPB PAGES/World Data 
Center-A for Paleoclimatology Data Contribution Series # 93-106. NOAA/NGDC Paleoclima- 
tology Program, Boulder, USA. 

Wessel P, Smith WHF (1998) New, improved version of generic mapping tools released. Eos. Trans. 
AGU 79:579. 

Wu P, Peltier WR (1983) Glacial isostatic adjustment and the free-air gravity anomaly as a 
constraint on deep mantle viscosity. Geophys. J. R. Astron. Soc. 74:377-449. 



37 



